This output report validates spatial pattern for the EU regions against LUH2v2 data.
# setup working dir and packages, read outputdir
readArgs("outputdir")
library(madrat)
library(mrcommons)
library(magpie4)
library(ggplot2)
library(plotly)
library(gridExtra)
library(patchwork)
library(dplyr)
library(tidyr)
library(stringr)
print(paste0("Script started for output directory: ", outputdir))
## [1] "Script started for output directory: ./output/v39k_FSECa_LandscapeElements_2024-05-24_16.53.30"
compareYears <- c(1995, 2000, 2005, 2010)
# ----- Read model output and luh2v2 data, process and bind them to a validation object -----
#### load and process reference data to match with model data
luh2v2 <- calcOutput("LanduseInitialisation", aggregate = FALSE, cellular = TRUE,
nclasses = "seven", cells = "lpjmlcell")
luh2v2 <- luh2v2[, compareYears, ]
# add subdimension for reference data
getNames(luh2v2) <- paste0(getNames(luh2v2), ".ref")
#### load and process model data to match with reference data
path2landoutput <- file.path(outputdir, "cell.land_0.5.mz")
cellLand <- read.magpie(path2landoutput)
cellLand <- cellLand[, compareYears, ]
# add subdimension for model data
getNames(cellLand) <- paste0(getNames(cellLand), ".mod")
load(file.path(outputdir, "spatial_header.rda"))
##### bind model and reference data
validationObj <- mbind(cellLand, luh2v2)
##### extract EU countries
mapping <- toolGetMapping("regionmappingH12.csv")
countriesEU <- mapping$CountryCode[mapping$RegionCode == "EUR"]
# ----- Create a plotting functions -----
convertToWideDataframe <- function(magclassObj) {
as.data.frame(magclassObj) %>%
pivot_wider(names_from = "Data2", values_from = "Value")
}
subsetCountryCluster <- function(magclassObj, countries) {
magclassObj[intersect(countries, getItems(magclassObj, split=TRUE)[[1]][[1]]), , ]
}
# combines a scatter plot with maps such that facets are aligned
plotCombined <- function(validationObj) {
# create the plots
p1 <- plotScatter(validationObj)
p2 <- plotMaps(validationObj)
# combine and align the plots using patchwork
combinedPlot <- p1 / p2
# print the combined plot
print(combinedPlot)
}
# create scatter plot
plotScatter <- function(validationObj) {
validationObj <- toolCoord2Isocell(validationObj, "lpjcell")
validationObj <- subsetCountryCluster(validationObj, countriesEU)
### Create quality measure dataframe with RMSE, MAE a.o.
qualityMeasuresDf <- NULL # list of dataframes, one for each year
numericYears <- as.numeric(gsub("y", "", getItems(validationObj, 2)))
for (year in numericYears) {
# get a vector with the quality measures
qualityMeasures <- luplot::qualityMeasure(
validationObj[, year, "mod"], validationObj[, year, "ref"],
measures = c("Willmott", "Willmott refined", "Nash Sutcliffe", "RMSE", "MAE")
)
# Create text for the year facet
text <- ""
for (i in seq_along(qualityMeasures)) {
text <- paste0(text, stringr::str_trunc(names(qualityMeasures)[i], 12, ellipsis = "."),
" : ", qualityMeasures[i], "\n")
}
# Create df for that year
qualityMeasuresDf[[year]] <- data.frame(
Year = year,
text = text
)
}
qualityMeasBinded <- dplyr::bind_rows(qualityMeasuresDf) # bind all year dataframes
# Create the plot with facets
validationDF <- convertToWideDataframe(validationObj)
p <- ggplot(validationDF,
aes(x = mod, y = ref, reg = Region, cell = Cell, relativeErr = (mod - ref) / ref)) +
geom_point(size = 0.5) +
geom_abline(color = "#663939", size = 1.5) +
facet_wrap(~ Year, scales = "free", ncol = 4) + # Create facets based on 'Year'
labs(x = "MAgPIE output", y = "luh2v2") +
theme(panel.background = element_rect(fill = "gray55"),
panel.grid.major = element_line(color = "gray62"),
panel.grid.minor = element_line(color = "grey58"))
# Add quality measure text
p <- p + geom_text(x = 0, y = Inf, aes(label = text), color = "#ffb5b5", # bright color contrasting dark background
hjust = 0, vjust = 1.1, nudge_x = 10, size = 2.6, family = "sans", fontface = "bold",
data = qualityMeasBinded, inherit.aes = FALSE)
return(p)
}
# create maps
plotMaps <- function(validationObj) {
# create a ggplot object using luplot
relErr <- (validationObj[, , paste0("mod")] - validationObj[, , paste0("ref")]) /
(validationObj[, , paste0("ref")])
p <- luplot::plotmap2(relErr,
legend_range = c(-2, 2), legendname = "relative\n diff \n to \n LUH2v2", ncol = 4,
midcol = "#ffffff", lowcol = "blue", highcol = "red", midpoint = 0,
title = ""
)
# adjust the plot
p <- p + coord_sf(xlim = c(-10, 40), ylim = c(35, 70)) + facet_wrap(~ Year, ncol = 4) +
theme(aspect.ratio = 1, legend.title = element_text(size = 8))
p <- p + guides(fill = guide_colorbar(barheight = 5, barwidth = 0.2))
return(p)
}
# create and interactive scatterplot for eu countries of a year
plotInteractiveSub <- function(validationObj, year) {
p <- plotScatter(validationObj[, year, ])
p <- p + ggtitle(paste0("EU Countries ", year))
p <- plotly::ggplotly(p)
p
}
type <- "crop"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)
type <- "past"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)
type <- "forestry"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)
type <- "primforest"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)
type <- "secdforest"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)
type <- "urban"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)
type <- "other"
plotCombined(validationObj[, , type])
plotInteractiveSub(validationObj[, , type], 2010)